######################################################################
################ ENVIRONMENT SETUP ###################################
######################################################################

# cleans the environment
rm(list=ls())

setwd("")

# loads the packages
library(corrplot)
library(vioplot)
library(reshape)
library(plotrix)
library(RColorBrewer)
library(plotly)
library(foreign)


colmap <- function(val, classint="quantile", palette="Blues", userbrk=NULL, nbrks=5, invpalet=F,custcol=NULL, transp=99, na.col=NULL){
  
  require(RColorBrewer)
  require(classInt)
  
  if(is.null(na.col)) val[is.na(val)] <- 0
  
  
  if (! is.null(custcol)){ 
    palette <- custcol
  }else{  
    pal <- palette
    palette <- brewer.pal(9, palette)
    if (pal=="Greys") palette <- palette[1:length(palette)-1]
    if (invpalet==T) palette <- palette[length(palette):1]
  }
  
  col <- colorRampPalette(palette)
  
  if (classint=="unclassed"){
    groups <- 0
    colors <- col(length(val))
    palette <- colors
    colors <- colors[findInterval(val,sort(val), all.inside=T)]
    pct.cases <- 0
  }
  
  
  if (classint=="unique"){
    groups
    groups <- sort(unique(val))
    
    if(is.character(val)) val <- as.factor(val)
    if(is.factor(val)) val <- as.numeric(val)
    if (length(groups)>length(custcol)) colors <- col(length(groups)) else colors <- custcol
    palette <- unique(colors)
    colors <- colors[findInterval(val,sort(as.numeric(unique(val))), all.inside=F)]
    pct.cases <- table(colors)[order(palette)]/length(colors)*100
  }
  
  
  
  if (classint=="user"){
    if (is.null(userbrk)) stop("El mapa de tipo 'User' requiere que los intervalos sean providos por el usuario!")
    
    groups <- userbrk
    colors <- col(length(userbrk)-1)
    palette <- unique(colors)
    colors <- colors[findInterval(val,userbrk, all.inside=T)] 
    pct.cases <- table(colors)[order(palette)]/length(colors)*100
    
  }
  
  if (classint %in% c("sd","quantile","equal","kmeans","jenks","pretty")){
    
    if (classint=="jenks"){
      if (length(val)>900){groups <- classIntervals(var=sample(val,120), n=nbrks, style=classint)
      }else {groups <- classIntervals(var=val, n=nbrks, style=classint)}
    }       
    if (classint!="jenks") groups <- classIntervals(var=val, n=nbrks, style=classint)
    
    groups <- unique(groups[[2]])
    colors <- col(length(groups)-1)
    palette <- unique(colors)
    colors <- colors[findInterval(val,groups, all.inside=T)] 
    pct.cases <- table(colors)[order(palette)]/length(colors)*100
    
  }
  
  if(length(colors[is.na(colors)])>0) colors[is.na(colors)] <- na.col
  
  
  if (transp!=99){
    colors <- paste0(colors, transp, sep="")
    palette <- paste0(palette, transp)}
  
  return(list(colors=colors, groups=groups, palette=palette, pct.cases=pct.cases))
  
}  



adleg <- function(legpos="bottomleft",leg.main="", cex.leg=1.4, pt.cex=0.8, y.intersp=0.8, x.intersp=0.5, pch=19, classint="quantile", val=NULL, palette=NULL,userbrk=NULL,groups=NULL,dec=1, plot.pct=FALSE, pct.cases=NULL, cl="SpatialPolygonsDataFrame",   big.mark=".", decimal.mark=",",horiz=F, add=T, symbol.bg="red", symbol.fg="black", leg.col="black", na.col=NULL, na.lab="n.a.", leg.bord="transparent"){
  
  require(maptools)
  
  par(mar=rep(0,4))
  
  if (is.numeric(legpos)) legpos <- cbind(legpos[1], legpos[2])
  
  if (classint=="unclassed"){
    
    # Open the library fields required to the progressive color bar
    require(fields)
    limites <- c(min(val,na.rm=T),max(val, na.rm=T))
    image.plot(legend.only=T, zlim=limites, col= sort(unique(palette),decreasing=T), smallplot=c(0.05,0.08,0.05,0.5))
  }  
  
  
  if (classint!="unclassed"){
    
    if (classint=="user"& is.null(userbrk)) stop("O mapa de tipo 'user' exige que os intervalos sejam fornecidos pelo usuário!")
    
    if (classint=="unique"){
      if (is.factor(val)) legs <- levels(val) else legs <- as.character(groups)}
    
    if (classint !="unique") legs <- leglabs(format(round(groups, dec), big.mark=big.mark, decimal.mark=decimal.mark), under="<", over=">")
    
    if (plot.pct==TRUE)legs <- paste(legs,paste("(", format(round(pct.cases, dec), big.mark=big.mark, decimal.mark=decimal.mark), "%)", sep=""), sep=" ")
    
    if(add==F) plot(1, col="transparent", axes=F)
    
    if (! is.null(na.col)) {
      legs <-  c(na.lab,legs)
      palette <- c(na.col, palette) 
    }
    
    if (cl%in% c("SpatialPolygons","SpatialPolygonsDataFrame")){
      legend(legpos,legend=legs,fill=palette, bg="transparent", bty="n", cex=cex.leg, y.intersp=y.intersp, x.intersp=x.intersp, title=leg.main, title.adj=0.15, horiz = horiz, text.col = leg.col, border = leg.bord)}
    else if (cl%in% c("SpatialPoints","SpatialPointsDataFrame")){
      legend(legpos,legend=legs, pt.bg = symbol.bg, border=symbol.fg, bg="transparent", bty="n", pch=pch, pt.cex=pt.cex, cex=cex.leg, y.intersp=y.intersp, x.intersp=x.intersp, title=leg.main, title.adj=0.15, horiz = horiz, text.col = leg.col, col=palette)
    }
    
  }
  
}

colBright <- function(colors, lim=130, bright="#ffffff", dark="#000000"){
  
  require(grDevices)
  
  rgb <- col2rgb(colors)
  
  for(i in 1:ncol(rgb)){
    val <- 0.2126*rgb[1,i]+0.7152*rgb[2,i]+0.0722*rgb[3,i]
    if(i==1) valc <- val else valc <- c(valc,val)
  }
  
  col <- rep(bright, length(colors))
  col[valc>lim] <- dark
  
  col
  
}

adlab <- function(spobj, labels, cex=1, font="Times New Roman", text.col="black", pos=NULL, col.bright=NULL, lab.font=1){
  
  require(sp)
  
  windowsFonts(times=windowsFont(font))
  par(family="times")
  
  if (! is.null(col.bright)){
    
    text.col <- colBright(col.bright)
    
  }
  
  if (! class(spobj) %in% c("SpatialPoints","SpatialPointsDataFrame")){
    text(coordinates(spobj), labels=labels, pos = pos, cex = cex, col = text.col, adj = 0.5, offset = 0,font=lab.font)
  }else{ 
    text(spobj, labels=labels, pos = pos, cex = cex, col = text.col, adj = 0.5,offset = 0,font=lab.font)
  }
  
}


mapa <- function(spobj=NULL, val=NULL, classint="quantile", nbrks=5, userbrk=NULL, title=NULL, palette="Blues", invpalet=F, legpos="bottomleft", dec=1,  cex.leg=1.4, y.intersp=0.5, x.intersp=0.5, smallplot=c(0.05,0.08,0.05,0.5),  border="black", lwd0=1, frame1=NULL, bordframe1="black",lwd1=1,frame2=NULL, bordframe2="black", lwd2=2, frame3=NULL, bordframe3="black", lwd3=2,pch=19, pt.cex=0.8, leg.main="", pct.cases=FALSE, custcol=NULL, transp=99, font="Times New Roman",leg.return=F, leg.draw=T, leg.horiz=F, big.mark=".", decimal.mark=",", labels=NULL,lab.cex=1, lab.col="black", lab.pos=NULL, lab.font=1, col.bright=NULL, add=F, leg.col="black",leg.bord="transparent",na.col=NULL,na.lab="n.a.",xlim=NULL, ylim=NULL,...){
  

  #   # Parte I - configuracao basica e checagem geral
  if (is.null(spobj) | is.null(val)) stop("Un valor válido debe ser suministrado para los parámetros mapobj (objeto espacial) y val (variable que se desea representar en el mapa)!")
  #   
  if(is.null(title)) par(mar=rep(0.1,4)) else par(mar=c(0.1,0.1,2,0.1))
  # if(! is.factor(val)) val[is.na(val)] <- 0
  
  require(maptools)
  
  if(class(val) %in% c("character","factor")){
    classint <- "unique"
    if (palette=="Blues") palette <- "Set3"}  
  
  
  
  # Parte II - obtem grupos e cores  
  mapres <- colmap(val = val, classint=classint, nbrks=nbrks, palette=palette, userbrk=userbrk, invpalet=invpalet, custcol=custcol, transp = transp, na.col = na.col)
  colores <- mapres[[1]]
  grupos <- mapres[[2]]
  paleta <- mapres[[3]] 
  pct.casos <- mapres[[4]]
  
  
  
  if (! is.null(frame1)){
    xlia <- frame1@bbox[1,]
    ylia <- frame1@bbox[2,]
  }else{
    xlia <- spobj@bbox[1,]
    ylia <- spobj@bbox[2,]
  }
  
  if(! is.null(xlim)) xlia <- xlim
  if(! is.null(ylim)) ylia <- ylim
  
  
  cl <- class(spobj)
  
  # Parte III - Desenha o mapa (e até três camadas superiores)
  if (cl %in% c("SpatialPolygons","SpatialPolygonsDataFrame")){
    plot(spobj, col=colores, border=border, lwd=lwd0, add=add,xlim=xlia, ylim=ylia)
  }
  else if (cl %in% c("SpatialPoints","SpatialPointsDataFrame")){
    plot(spobj, col=colores,cex=pt.cex,pch=pch, add=add,xlim=xlia, ylim=ylia)
  }
  
  if (! is.null(frame1)) plot(frame1, border=bordframe1, add=T, lwd=lwd1,xlim=xlia, ylim=ylia)
  if (! is.null(frame2)) plot(frame2, border=bordframe2, add=T, lwd=lwd2,xlim=xlia, ylim=ylia)
  if (! is.null(frame3)) plot(frame3, border=bordframe3, add=T, lwd=lwd3,xlim=xlia, ylim=ylia)
  
  #   # Parte IV: Adiciona o t??tulo (se existente), legenda ou barra de cores (se "unclassed")
  if (! is.null(title)) title(main=title)
  
  
  if (! is.null(labels)) adlab(spobj=spobj, labels=labels, cex=lab.cex, font=font, text.col=lab.col, pos=lab.pos, col.bright=colores,lab.font = lab.font)
  
  if (leg.return==T) return(list(legpos=legpos,leg.main=leg.main, cex.leg=cex.leg, y.intersp=y.intersp,x.intersp=x.intersp, pch=pch, classint=classint, val=val, palette=paleta,userbrk=userbrk,groups=grupos,dec=dec, pct.cases=pct.cases, pct.casos=pct.casos, cl=cl, leg.col=leg.col))
  
  if (leg.draw==T) adleg(legpos = legpos,leg.main = leg.main, cex.leg = cex.leg, y.intersp = y.intersp,x.intersp = x.intersp, pch = pch, classint = classint, val = val, palette = paleta,userbrk = userbrk,groups = grupos,dec=dec, plot.pct = pct.cases, pct.cases = pct.casos, cl = cl, big.mark = big.mark, decimal.mark = decimal.mark, horiz=leg.horiz, leg.col=leg.col, na.col=na.col, na.lab=na.lab, leg.bord=leg.bord)
  
}



mvioplot <- function (x, by=NULL, range = 1.5, h = NULL, ylim = NULL, names = NULL, horizontal = FALSE, col = "black", lty = 1, lwd = 1, rectCol = "gray50", rectBord = "gray50", colMed = "black", colMean="red", pchMed = 19, pchMean=19, at, add = FALSE, wex = 1, drawRect = TRUE, side="both", box=F, lag=0.03, plot.median=T, plot.mean=F, lwd.line = 1, lty.line = 1, col.line="black", inv.y=F, las=1){
  
  require(sm)
  
  side <- tolower(side)
  
  if (is.null(names) & ! is.null(by)) names <- unique(as.character(by))
  
  if (side%in%c("up","down","hboth")) horizontal <- TRUE
  
  if(! is.null(by)) datas <- split(x,by) else datas <- x
  
  
  
  if(inv.y==T) datas <- datas[sort(names(datas),decreasing = T)]
  
  if (class(datas)== "list") n <- length(datas) else n <- 1
  
  if (missing(at)){
    at <- 1:n
    
    if (side%in%c("left","down")) at <- at - lag else at <- at + lag 
    
  }  
  
  upper <- vector(mode = "numeric", length = n)
  lower <- vector(mode = "numeric", length = n)
  q1 <- vector(mode = "numeric", length = n)
  q3 <- vector(mode = "numeric", length = n)
  med <- vector(mode = "numeric", length = n)
  mea <- vector(mode = "numeric", length = n)
  base <- vector(mode = "list", length = n)
  height <- vector(mode = "list", length = n)
  baserange <- c(Inf, -Inf)
  args <- list(display = "none")
  
  if (!(is.null(h))) args <- c(args, h = h)
  
  for (i in 1:n) {
    if(class(datas)=="list") data <- datas[[i]] else data <- datas
    data <- data[! is.na(data) & ! is.infinite(data) & ! is.nan(data)]
    data.min <- min(data, na.rm = T)
    data.max <- max(data, na.rm = T)
    q1[i] <- quantile(data, 0.25, na.rm = T)
    q3[i] <- quantile(data, 0.75, na.rm = T)
    med[i] <- median(data, na.rm = T)
    mea[i] <- mean(data, na.rm = T)
    iqd <- q3[i] - q1[i]
    upper[i] <- min(q3[i] + range * iqd, data.max, na.rm = T)
    lower[i] <- max(q1[i] - range * iqd, data.min, na.rm = T)
    est.xlim <- c(min(lower[i], data.min,na.rm = T), max(upper[i], data.max, na.rm = T))
    smout <- do.call("sm.density", c(list(data, xlim = est.xlim), args))
    
    
    hscale <- 0.4/max(smout$estimate) * wex
    base[[i]] <- smout$eval.points
    height[[i]] <- smout$estimate * hscale
    t <- range(base[[i]])
    baserange[1] <- min(baserange[1], t[1])
    baserange[2] <- max(baserange[2], t[2])
  }
  
  
  if (!add){
    xlim <- if (n == 1) at + c(-0.5, 0.5) else range(at) + min(diff(at))/2 * c(-1, 1)
    if (is.null(ylim)) {
      ylim <- baserange
      
    }
  }
  if (is.null(names)) {
    label <- 1:n
  }else{
    label <- names
  }
  boxwidth <- 0.05 * wex
  
  if (!add) plot.new()
  
  if (!horizontal) {
    if (!add) {
      plot.window(xlim = xlim, ylim = ylim)
      axis(2, las=las)
      axis(1, at = at, label = label)
    }
    
    if (box==T) box()
    
    for (i in 1:n) {
      
      
      if(side%in%c("left","down")){
        cvl <- c(at[i] - height[[i]])
        lines(x=cvl, y=base[[i]], col = col, lty = lty, lwd = lwd)
      }else if (side%in%c("right","up")){
        cvr <- c(at[i] + height[[i]])
        lines(x=rev(cvr), y=rev(base[[i]]), col = col, lty = lty, lwd = lwd)
      }else{
        cvl <- c(at[i] - height[[i]])
        cvr <- c(at[i] + height[[i]])
        lines(x=cvl, y=base[[i]], col = col, lty = lty, lwd = lwd)
        lines(x=rev(cvr), y=rev(base[[i]]), col = col, lty = lty, lwd = lwd)
      }
      
      if (drawRect) {
        lines(at[c(i, i)], c(lower[i], upper[i]), lwd = lwd.line, lty = lty.line, col=col.line)
        rect(at[i] - boxwidth/2, q1[i], at[i] + boxwidth/2, q3[i], col = rectCol, border = rectBord)
        
        if (plot.median==T) points(at[i], med[i], pch = pchMed, col = colMed)
        
        if (plot.mean==T) points(at[i], mea[i], pch = pchMean, col = colMean)
        
      }
    }
  }else{
    if (!add){
      plot.window(xlim = ylim, ylim = xlim)
      axis(1)
      axis(2, at = at, label = label, las=las)
    }
    
    if (box==T) box()    
    
    for (i in 1:n){ 
      
      if(side%in%c("left","down")){
        cvl <- c(at[i] - height[[i]])
        lines(x=base[[i]], y=cvl, col = col, lty = lty, lwd = lwd)
      }else if (side%in%c("right","up")){
        cvr <- c(at[i] + height[[i]])
        lines(x=rev(base[[i]]), y=rev(cvr), col = col, lty = lty, lwd = lwd)
      }else{
        cvl <- c(at[i] - height[[i]])
        cvr <- c(at[i] + height[[i]])
        lines(x=base[[i]], y=cvl, col = col, lty = lty, lwd = lwd)
        lines(x=rev(base[[i]]), y=rev(cvr), col = col, lty = lty, lwd = lwd)
      }
      
      
      # polygon(c(base[[i]], rev(base[[i]])), c(at[i] - height[[i]], 
      #                                         rev(at[i] + height[[i]])), col = col, border = border, 
      #         lty = lty, lwd = lwd)
      if (drawRect) {
        lines(c(lower[i], upper[i]), at[c(i, i)], lwd = lwd.line, lty = lty.line, col=col.line)
        rect(q1[i], at[i] - boxwidth/2, q3[i], at[i] + boxwidth/2, col = rectCol, border = rectBord)
        
        if (plot.median==T) points(med[i], at[i], pch = pchMed, col = colMed)
        
        if (plot.mean==T) points(mea[i], at[i], pch = pchMean, col = colMean)
      }
    }
  }
  invisible(list(upper = upper, lower = lower, median = med, 
                 q1 = q1, q3 = q3))
}


geraSample <- function (var, w, by, data){
  
  byv <- unique(data[,by])
  
  for(i in 1:length(byv)){
    
    dt <- data[data[,by]==byv[i],]
    
    sp <- sample(rep(dt[,var],dt[,w]), size = length(dt[,var]), replace = T)
    sp <- cbind(byv[i],sp)  
    
    if(i==1) st <- sp else st <- rbind(st,sp)
    
  }
  
  st <- data.frame(st)
  names(st) <- c("by","var")    
  
  st  
  
}  

dfDesc <- function(var,dec=1, big.mark="", return=F){
  require(pastecs)
  
  desc <- stat.desc(var)
  
  n <- format(round(desc[1]), big.mark = ",")
  nas <- format(round(desc[3]), big.mark = ",")
  mini <- format(round(desc[4], dec), big.mark = big.mark, small.mark = ".")
  maxi <- format(round(desc[5], dec), big.mark = big.mark, small.mark = ".")
  rang <- format(round(desc[6], dec), big.mark = big.mark, small.mark = ".")
  med <- format(round(desc[8], dec), big.mark = big.mark, small.mark = ".")
  mea <- format(round(desc[9], dec), big.mark = big.mark, small.mark = ".")
  ste <- format(round(desc[13], dec), big.mark = big.mark, small.mark = ".")
  cfv <- format(round(desc[14], dec), big.mark = big.mark, small.mark = ".")
  
  dt <- data.frame(N=n,NAs=nas,Min=mini,Max=maxi,Range=rang,Median=med,Mean=mea,Std.Dev.=ste,Coef.Var=cfv)
  row.names(dt)<- ""
  
  if(return==T) return(dt) else knitr::kable(dt)  
  
}


splag <- function(spobj, val, queen=T){
  require(spatstat)
  require(spdep)
  vecino <- poly2nb(spobj, queen=queen)
  lwvec <- nb2listw(vecino, zero.policy=T, )
  vallag <- lag.listw(lwvec, val, NAOK=T)
  vallag[is.na(vallag)] <- val[is.na(vallag)] 
  return(vallag)
}


plotConf <- function(x,y,error=0,line.col="black",add=F, opacity=10, ylim=c(0,40), nsp=50){
  
  require(gplots)
  
  ly <- y-error
  hy <- y+error
  
  sy <- spline(y, n = nsp)
  sly <- spline(ly, n = nsp)
  shy <- spline(hy, n = nsp)
  
  
  if (add==F){
    # par(mar=c(4,6,1,1))
    plot(sy, type="l", bty="n", ylim=ylim, xlab="", ylab="", axes = F)
    axis(1, at = c(1:length(x)), labels = x, col.ticks = "grey90", col = "grey90")
    # axis(2, at = seq(0,40,10), labels = seq(0,40,10), col.ticks = "grey90", col = "grey90", las=1)
    axis(2, at = seq(-0.5,0.5,0.5), labels = c("Center-Left","Center","Center-Right"), col.ticks = "grey90", col = "grey90", las=1)
    
    polygon(c(sy$x,rev(sy$x)),c(sly$y,rev(shy$y)),col = paste0(col2hex(line.col),opacity), border = FALSE)
    lines(sy, lwd = 2, col=line.col)
    
  }else{
    polygon(c(sy$x,rev(sy$x)),c(sly$y,rev(shy$y)),col = paste0(col2hex(line.col),opacity), border = FALSE)
    lines(sy, lwd = 2, col=line.col)
  }
  
}

is.even <- function(x) x %% 2 == 0 

tfscale <- function(x){
  round((x-min(x,na.rm = T))/(max(x,na.rm = T) - min(x,na.rm = T)) * 100,1)  
  
}


# loads the data
ef <- read.table("Data/BPSR_Data_Power_Rodrigues-Silveira_v2.csv", sep=";", dec=",", header=T)

######################################################################
########### Descriptives used in Fig.1 ###############################
######################################################################


a <- dfDesc(ef$ideo_imp, dec = 3, return = T)
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==1994], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==1996], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==1998], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2000], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2002], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2004], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2006], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2008], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2010], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2012], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2014], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2016], dec = 3, return = T))
a <- rbind(a,dfDesc(ef$ideo_imp[ef$year==2018], dec = 3, return = T))

a$year <- c("All years", seq(1994,2018,2))

a <- a[,c("year","Mean","Median","Std.Dev.","Min","Max")]
names(a) <- c("year","mean","median","sd","min","max")

tb <- a

rm(a)

for(i in 2:ncol(tb)) tb[,i] <- as.numeric(as.character(tb[,i]))

tb

######################################################################
########### Figure 1 - Margin Plot ###################################
######################################################################

par(mar=c(4,6,0,1))
plotConf(tb$year, tb$mean, error = tb$sd, ylim = c(-0.1,0.8), nsp = 50, add=F)
abline(h=0, lwd=1, col="red", lty=2)



######################################################################
########### Figure 2. MAP PANNEL #####################################
######################################################################


library(maptools)



load("Data/Municipios_SHP.RData")

year  <- seq(1994,2018,2)
inc <- c("Itamar","FHC","FHC","FHC","FHC","Lula","Lula","Lula","Lula","Dilma","Dilma","Temer","Temer")


png(filename = "Images/Map_Panel_v3.png", width = 3912, height = 4480)

par(mfrow=c(5,3))

for (i in 1:length(year)){
  
  af <- ef[ef$year==year[i],]
  af <- unique(af)
  
  if(year[i]%in%c(1994:2002)) b94 <- merge(b00, af, by="GEOCODIG_M", all.x=T)
  if(year[i]==2004) b94 <- merge(b05, af, by="GEOCODIG_M", all.x=T)
  if(year[i]%in%c(2006:2010)) b94 <- merge(b10, af, by="GEOCODIG_M", all.x=T, duplicateGeoms=T)
  if(year[i]>2010) b94 <- merge(b15, af, by="GEOCODIG_M", all.x=T)
  
  retleg <- mapa(b94, splag(b94,b94$ideo_na), border="transparent", frame2 = u, dec=3, decimal.mark = ".", big.mark = ",", y.intersp = 0.7, cex.leg = 1, classint = "user", userbrk = c(-1,0,0.15,0.25,0.40,1) , palette = "Greys", na.col = "gray", invpalet = F, xlim=c(-75.55451,-31.10946), ylim=c(-35.995777,7.235437), leg.draw=F, leg.return=T,  lwd2=4)
  
  
  if (is.even(i)) text(-40,2, labels = year[i], cex = 12, font = 2) else text(-40,2, labels = year[i], cex = 12)
  
  mu <- length(b94$ideo_na[b94$ideo_na<0])/length(b94$ideo_na)
  ri <- (( (mu*5)/5)^0.5)*5

  symbols(x=-66.6,y=-28.9+5,circle=5,inches=FALSE,add=TRUE)
  symbols(x=-66.6,y=-28.9+ri,circle=ri,inches=FALSE,add=TRUE, bg = "#F7F7F7")
  
  text(x=-66.6,y=-28.9+(ri*2)+2, labels = paste0(round(mu*100,1),"%"), cex=6)
  
}

par(mfrow=c(1,1))

dev.off()


######################################################################
########### Figure 3. President, CD, MIS #############################
######################################################################


a <- cbind("All",mean(ef$ideo_imp, na.rm=T),weighted.mean(ef$ideo_imp, ef$qtde_eleitores, na.rm=T))
a <- rbind(a,cbind(1994,mean(ef$ideo_imp[ef$year==1994]), weighted.mean(ef$ideo_imp[ef$year==1994], ef$qtde_eleitores[ef$year==1994])))
a <- rbind(a,cbind(1996,mean(ef$ideo_imp[ef$year==1996]), weighted.mean(ef$ideo_imp[ef$year==1996], ef$qtde_eleitores[ef$year==1996])))
a <- rbind(a,cbind(1998,mean(ef$ideo_imp[ef$year==1998]),weighted.mean(ef$ideo_imp[ef$year==1998], ef$qtde_eleitores[ef$year==1998])))
a <- rbind(a,cbind(2000,mean(ef$ideo_imp[ef$year==2000]),weighted.mean(ef$ideo_imp[ef$year==2000], ef$qtde_eleitores[ef$year==2000])))
a <- rbind(a,cbind(2002,mean(ef$ideo_imp[ef$year==2002], na.rm = T),weighted.mean(ef$ideo_imp[ef$year==2002], ef$qtde_eleitores[ef$year==2002], na.rm = T)))
a <- rbind(a,cbind(2004,mean(ef$ideo_imp[ef$year==2004], na.rm = T),weighted.mean(ef$ideo_imp[ef$year==2004], ef$qtde_eleitores[ef$year==2004], na.rm = T)))
a <- rbind(a,cbind(2006,mean(ef$ideo_imp[ef$year==2006], na.rm = T),weighted.mean(ef$ideo_imp[ef$year==2006], ef$qtde_eleitores[ef$year==2006], na.rm = T)))
a <- rbind(a,cbind(2008,mean(ef$ideo_imp[ef$year==2008]),weighted.mean(ef$ideo_imp[ef$year==2008], ef$qtde_eleitores[ef$year==2008])))
a <- rbind(a,cbind(2010,mean(ef$ideo_imp[ef$year==2010]),weighted.mean(ef$ideo_imp[ef$year==2010], ef$qtde_eleitores[ef$year==2010])))
a <- rbind(a,cbind(2012,mean(ef$ideo_imp[ef$year==2012], na.rm = T),weighted.mean(ef$ideo_imp[ef$year==2012], ef$qtde_eleitores[ef$year==2012], na.rm = T)))
a <- rbind(a,cbind(2014,mean(ef$ideo_imp[ef$year==2014]),weighted.mean(ef$ideo_imp[ef$year==2014], ef$qtde_eleitores[ef$year==2014])))
a <- rbind(a,cbind(2016,mean(ef$ideo_imp[ef$year==2016], na.rm = T),weighted.mean(ef$ideo_imp[ef$year==2016], ef$qtde_eleitores[ef$year==2016])))
a <- rbind(a,cbind(2018,mean(ef$ideo_imp[ef$year==2018], na.rm = T),weighted.mean(ef$ideo_imp[ef$year==2018], ef$qtde_eleitores[ef$year==2018])))


b <- cbind("all",mean(ef$ideo_na),weighted.mean(ef$ideo_na, ef$qtde_eleitores))

b <- rbind(b,cbind(1994,mean(ef$ideo_na[ef$year==1994]),weighted.mean(ef$ideo_na[ef$year==1994], ef$qtde_eleitores[ef$year==1994])))

b <- rbind(b,cbind(1996,mean(ef$ideo_na[ef$year==1996]), weighted.mean(ef$ideo_na[ef$year==1996], ef$qtde_eleitores[ef$year==1996])))

b <- rbind(b,cbind(1998,mean(ef$ideo_na[ef$year==1998]),weighted.mean(ef$ideo_na[ef$year==1998], ef$qtde_eleitores[ef$year==1998])))

b <- rbind(b,cbind(2000,mean(ef$ideo_na[ef$year==2000]),weighted.mean(ef$ideo_na[ef$year==2000], ef$qtde_eleitores[ef$year==2000])))

b <- rbind(b,cbind(2002,mean(ef$ideo_na[ef$year==2002], na.rm = T),weighted.mean(ef$ideo_na[ef$year==2002], ef$qtde_eleitores[ef$year==2002], na.rm = T)))

b <- rbind(b,cbind(2004,mean(ef$ideo_na[ef$year==2004], na.rm = T),weighted.mean(ef$ideo_na[ef$year==2004], ef$qtde_eleitores[ef$year==2004], na.rm = T)))

b <- rbind(b,cbind(2006,mean(ef$ideo_na[ef$year==2006]),weighted.mean(ef$ideo_na[ef$year==2006], ef$qtde_eleitores[ef$year==2006])))

b <- rbind(b,cbind(2008,mean(ef$ideo_na[ef$year==2008]),weighted.mean(ef$ideo_na[ef$year==2008], ef$qtde_eleitores[ef$year==2008])))

b <- rbind(b,cbind(2010,mean(ef$ideo_na[ef$year==2010]),weighted.mean(ef$ideo_na[ef$year==2010], ef$qtde_eleitores[ef$year==2010])))

b <- rbind(b,cbind(2012,mean(ef$ideo_na[ef$year==2012], na.rm = T),weighted.mean(ef$ideo_na[ef$year==2012], ef$qtde_eleitores[ef$year==2012], na.rm = T)))

b <- rbind(b,cbind(2014,mean(ef$ideo_na[ef$year==2014]),weighted.mean(ef$ideo_na[ef$year==2014], ef$qtde_eleitores[ef$year==2014])))

b <- rbind(b,cbind(2016,mean(ef$ideo_na[ef$year==2016], na.rm = T),weighted.mean(ef$ideo_na[ef$year==2016], ef$qtde_eleitores[ef$year==2016])))

b <- rbind(b,cbind(2018,mean(ef$ideo_na[ef$year==2018], na.rm = T),weighted.mean(ef$ideo_na[ef$year==2018], ef$qtde_eleitores[ef$year==2018])))


a <- data.frame(a)
names(a) <- c("year","mean_imp","w.mean_imp")
b <- data.frame(b)
names(b) <- c("year","mean_bls","w.mean_bls")

c <- merge(b, a, by="year" )
c$mean_bls <- round(as.numeric(as.character(c$mean_bls)),3)
c$w.mean_bls <- round(as.numeric(as.character(c$w.mean_bls)),3)
c$mean_imp <- round(as.numeric(as.character(c$mean_imp)),3)
c$w.mean_imp <- round(as.numeric(as.character(c$w.mean_imp)),3)

# c <- merge(c,ca, by="year")



# data with BLS values for the MIS
aga <- aggregate(ef$ideo_na, by=list(year=ef$year), mean, na.rm=T)
aga$ideo_wavg <- c$w.mean_bls[1:13]
aga$ideopres <- c(0.1935945255,0.2801891265,0.2801891265,0.2801891265,0.2801891265,-0.419610896,-0.419610896,-0.419610896,-0.419610896,-0.4231020425,-0.4231020425,0.570742088,0.570742088)
aga$popular <- c(0.5957,0.6125,0.61,0.4004,0.4740,0.6057,0.6286,0.7984,0.87,0.7614,0.5495,0.3950,0.12)
names(aga) <- c("year","ideo_avg","ideo_wavg","ideopres","popular")



# data with IMPUTED values for the MIS
agi <- aggregate(ef$ideo_imp, by=list(year=ef$year), mean, na.rm=T)
agi$ideo_wavg <- c$w.mean_imp[1:13]
agi$ideopres <- c(0.1935945255,0.2801891265,0.2801891265,0.2801891265,0.2801891265,-0.419610896,-0.419610896,-0.419610896,-0.419610896,-0.4231020425,-0.4231020425,0.570742088,0.570742088)
agi$popular <- c(0.5957,0.6125,0.61,0.4004,0.4740,0.6057,0.6286,0.7984,0.87,0.7614,0.5495,0.3950, 0.12)
names(agi) <- c("year","ideo_avg","ideo_wavg","ideopres","popular")

ag <- agi

ag$size <- (ag$popular/max(ag$popular))^0.57*20


se <- read.table("Data/party_seats.csv", header=T, sep=";")
id <- read.csv("Data/partidos_ano_v2018.csv", header=T, sep=";")

se <- merge(se, id, by.x=c("year","party"), by.y=c("ANO_ELEICAO","NUM_PART"), all.x=T)
se <- se[se$seats>0,]

se$prop_seats <- se$seats/513

se$ct_congress <- se$prop_seats*se$IDEOLOGY

se <- se[! is.na(se$IDEOLOGY),]

ap <-aggregate(se$ct_congress, by=list(year=se$year), sum)

ap$party <- 99

ap$accronym <- "CD"

names(ap) <- c("year","IDEOLOGY","party","accronym")
ap <- ap[,c(1,3,4,2)]

sp <- se[se$party%in%c(13,45,15,25), c("year","party","accronym","IDEOLOGY")]

sp <- rbind(sp,ap)

sp$accronym <- gsub(pattern = ".",replacement = "/", fixed = T, x = sp$accronym)


agw <- ag[ag$year%in% seq(1994,2018,4), c("year","ideo_wavg")]

agm <- ag[ag$year%in% seq(1994,2018,4), c("year","ideo_avg")]

pres_m <- c(-0.221789,-0.1378,-0.3899,-0.222,-0.2594,-0.2317,0.0484)
pres_wm <- c(-0.2375849, -0.2174, -0.4611,-0.2262,-0.242,-0.21536, 0.1682)


par(mfrow=c(2,2))
{
  # PLOT 1 - MIS
  
  par(mar=c(4,8,1,1))
  # MIS (Unweighted)
  plot(spline(agm$year, round(agm$ideo_avg,2), n=500), type="l", bty="n", axes=F, xlab="", ylab="", ylim=c(0,0.5), lwd=2, lty=2, col="grey80") 
  
  axis(1, at = seq(1994,2018,4),labels = seq(1994,2018,4))
  
  axis(2, at = c(0,0.5),labels = c("Center (0.0)","Center-Right (0.5)"), las=1)
  
  
  # MIS (Weighted)
  lines(spline(agw$year, round(agw$ideo_wavg,2), n=500), lty=1, col="black", lwd=2) 
  
  abline(h=0, col="grey50", lwd=2)
  
  legend(x = 2002, y=0.48, legend = c("MIS (Unweighted)","MIS (Weighted)"), bty="n", lty = c(2,1), lwd=2, col=c("grey80","black"))
  
  
  # PLOT 2 - Entrant President and CD
  
  par(mar=c(4,8,1,1))
  # President
  plot(spline(agm$year, round(pres_wm,2), n=500), type="l", bty="n", axes=F, xlab="", ylab="", ylim=c(-0.5,0.5), lwd=2, lty=1, col="black") 
  
  axis(1, at = seq(1994,2018,4),labels = seq(1994,2018,4))
  
  axis(2, at = c(-0.5,0,0.5),labels = c("Center-Left (-0.5)","Center (0.0)","Center-Right (0.5)"), las=1)
  
  
  # CD
  lines(spline(sp$year[sp$party==99], round(sp$IDEOLOGY[sp$party==99],2), n=500), lty=2, col="grey80", lwd=2) 
  
  abline(h=0, col="grey50", lwd=2)
  
  legend(x = 2002, y=0.5, legend = c("Chamber of Deputies","MIS (President)"), bty="n", lty = c(2,1), lwd=2, col=c("grey80","black"))
  
  
  
  # PLOT 3 - PT and PSDB
  
  par(mar=c(4,8,1,1))
  # PT
  plot(spline(sp$year[sp$party==13], round(sp$IDEOLOGY[sp$party==13],2), n=500), type="l", bty="n", axes=F, xlab="", ylab="", ylim=c(-1,0.5), lwd=2, lty=2, col="grey80") 
  
  axis(1, at = seq(1994,2018,4),labels = seq(1994,2018,4))
  
  axis(2, at = c(-1,-0.5,0,0.5),labels = c("Left (-1.0)","Center-Left (-0.5)","Center (0.0)","Center-Right (0.5)"), las=1)
  
  
  # PSDB
  lines(spline(sp$year[sp$party==45], round(sp$IDEOLOGY[sp$party==45],2), n=500), lty=1, col="black", lwd=2) 
  
  abline(h=0, col="grey50", lwd=2)
  
  legend(x = 2008, y=-0.6, legend = c("PT","PSDB"), bty="n", lty = c(2,1), lwd=2, col=c("grey80","black"))
  
  
  
  # PLOT 4 - MDB and DEM
  
  par(mar=c(4,8,1,1))
  # PMDB
  plot(spline(sp$year[sp$party==15], round(sp$IDEOLOGY[sp$party==15],2), n=500), type="l", bty="n", axes=F, xlab="", ylab="", ylim=c(0,1), lwd=2, lty=2, col="grey80") 
  
  axis(1, at = seq(1994,2018,4),labels = seq(1994,2018,4))
  
  axis(2, at = c(0,0.5,1),labels = c("Center (0.0)","Center-Right (0.5)","Right (1.0)"), las=1)
  
  
  # PFL/DEM
  lines(spline(sp$year[sp$party==25], round(sp$IDEOLOGY[sp$party==25],2), n=500), lty=1, col="black", lwd=2) 
  
  abline(h=0, col="grey50", lwd=2)
  
  legend(x = 1994, y=0.5, legend = c("(P)MDB","PFL/DEM"), bty="n", lty = c(2,1), lwd=2, col=c("grey80","black"))
  
  
}
par(mfrow=c(1,1))




######################################################################
########### Figure 4. President idelogy, aproval #####################
######################################################################



ag$size <- (ag$popular/max(ag$popular))^0.57*5

ag <- ag[order(ag$year, decreasing = T),]

sp <- spline(x = ag$ideo_avg[nrow(ag):1], y = ag$year[nrow(ag):1])
ag <- ag[order(ag$year),]

sp <- spline(x=ag$ideo_avg, y=ag$year)
sp <- sp[order(sp$y),]


par(mar=c(2,4,1,1))
# MIS (Unweighted)
plot(x=ag$ideo_avg, y=ag$year, type="l", bty="n", axes=F, xlab="", ylab="", xlim=c(-1,1), ylim=c(2018,1994),lwd=2, lty=2, col="grey80") 
axis(1, at = c(-1,0,1), labels = c("Left","Center","Right"))
axis(2, at = seq(2018,1994,-2), labels = seq(2018,1994,-2), las=1)

lines(x=ag$ideo_wavg, ag$year, lwd=2, lty=1, col="black")
abline(v=0, col="grey50", lwd=2)

points(x=ag$ideopres, y=ag$year, pch=19, cex=ag$size, col="#58585860")

text(y=1993.9, x=0.5, labels = "Itamar")
text(y=1995.9, x=0.55, labels = "FHC")
text(y=1997.9, x=0.55, labels = "FHC")
text(y=1999.9, x=0.55, labels = "FHC")
text(y=2001.9, x=0.55, labels = "FHC")
text(y=2003.9, x=-0.68, labels = "Lula")
text(y=2005.9, x=-0.68, labels = "Lula")
text(y=2007.9, x=-0.68, labels = "Lula")
text(y=2009.9, x=-0.7, labels = "Lula")
text(y=2011.9, x=-0.72, labels = "Dilma")
text(y=2013.9, x=-0.72, labels = "Dilma")
text(y=2015.9, x=0.65, labels = "Temer", adj=0)
text(y=2017.9, x=0.63, labels = "Temer", adj = 0)

text(-0.99, 1994, labels = "Ideological Score", adj=0, font = 2)

points(-0.95, 1995, cex=2, col="#58585860", pch=19, border="#58585860")
text(-0.91, 1994.9, labels = "Presidential Placement", adj=0)

legend(-1.05,1995, legend = c("MIS (Unweighted)", "MIS (Weighted)"), col = c("grey80","black"), lty=c(2,1), x.intersp = 0.5, bty="n", seg.len = 0.7)

text(-0.99, 1998, labels = "Presidential Approval", adj=0, font = 2)

legend(-1.05,1998.4, legend = c("12%","40%","87%"), col = c("grey80","grey80","grey80"), pch=19, pt.cex = c(1.616506,3.21,5), bty="n", horiz=T, x.intersp = c(0.5,1,1.4))






######################################################################
########### Figure A1. Multivariate Violin Plot  #####################
######################################################################


mvioplot(ef$ideo_na , ef$year, side="up", lag=0.05, ylim=c(-1,1), lwd=2, inv.y = T, las=1)
mvioplot(ef$ideo_imp , ef$year, side="up", lag=0.05, ylim=c(-1,1), lwd=2, inv.y = T, lty=2, drawRect = F, add=T, col = "grey80")

wag <- geraSample(var = "ideo_na",w = "qtde_eleitores",by = "year", data = ef) 


mvioplot(wag$var, wag$by, side="down", lag=0.05, ylim=c(-1,1), add=T, lwd=2, lty=3)

wag <- geraSample(var = "ideo_imp",w = "qtde_eleitores",by = "year", data = ef) 

mvioplot(wag$var, wag$by, side="down", lag=0.05, ylim=c(-1,1), col="grey80", drawRect = F, add=T, lwd=2, lty=4)

abline(v=0)

legend(y=8,x=-1, legend=rep("",4), col = c("black","grey80","black","grey80"), lty=c(1,2,3,4), bty="n", seg.len = 0.7)



######################################################################################################
######################################################################################################
######################################################################################################
########################### Regressions ##############################################################
######################################################################################################
######################################################################################################
######################################################################################################



rm(list=ls())


gc()
# Opens the data
ef <- read.csv("Data/BPSR_Data_Power_Rodrigues-Silveira_v2.csv", header=T, sep=";",dec = ",")

ef$ln_voter <- log(ef$qtde_eleitores)
ef$ln_voter[is.infinite(ef$ln_voter)] <- NA


for(i in 1:ncol(ef)){
  
  ef[,i][is.nan(ef[,i])] <- NA
  ef[,i][is.infinite(ef[,i])] <- NA
  
  
}


ef <- ef[! is.na(ef$closeness),]

ef$ln_voter <- log(ef$qtde_eleitores)
ef$ln_voter[is.infinite(ef$ln_voter)] <- NA

names(ef)[16] <- "fragf"

ez <- scale(ef[c(6:ncol(ef))])
ez <- data.frame(ez)




######################################################################################################
######################################################################################################
######################################################################################################
########################### OLS Regressions with the BLS D.V. AND LAGGED D.V. ########################
######################################################################################################
######################################################################################################
######################################################################################################


# Regression 1: imcumbent spillover effect
reg1 <- lm(ef$ideo_na~ef$ideo_na_lag+ef$closeness+ef$close_ideopres+ef$coal_pres+ef$ideo_pres+ef$supp_pres+I(ef$supp_pres*ef$ideo_pres)+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg1)


# Regression 2: socioeconomic modernization
reg2 <- lm(ef$ideo_na~ef$ideo_na_lag+ef$IDHM+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg2)


# Regression 3: political modernization (pluralism)
reg3 <- lm(ef$ideo_na~ef$ideo_na_lag+ef$fragf+ef$comp+ef$turnout+ef$pol_pi+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg3)


# Regression 4: pro-poor policies and PT pro-poor policies
reg4 <- lm(ef$ideo_na~ef$ideo_na_lag+ef$var_pob+ef$pt+I(ef$var_pob*ef$pt)+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg4)


# Regression 5: pooled regression
reg5 <- lm(ef$ideo_na~ef$ideo_na_lag+ef$closeness+ef$ideo_pres+ef$close_ideopres+ef$supp_pres+I(ef$supp_pres*ef$ideo_pres)+ef$coal_pres+ef$IDHM+ef$fragf+ef$comp+ef$turnout+ef$pol_pi+ef$var_pob+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg5)



######################################################################################################
######################################################################################################
######################################################################################################
########################### OLS Regressions with the IMPUTED D.V. AND LAGGED D.V. ####################
######################################################################################################
######################################################################################################
######################################################################################################



# Regression 1: imcumbent spillover effect
reg1i <- lm(ef$ideo_imp~ef$ideo_imp_lag+ef$closeness+ef$close_ideopres+ef$coal_pres+ef$ideo_pres+ef$supp_pres+I(ef$supp_pres*ef$ideo_pres)+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg1i)

# Regression 2: socioeconomic modernization
reg2i <- lm(ef$ideo_imp~ef$ideo_imp_lag+ef$IDHM+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg2i)


# Regression 3: political modernization (pluralism)
reg3i <- lm(ef$ideo_imp~ef$ideo_imp_lag+ef$fragf+ef$comp+ef$turnout+ef$pol_pi+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg3i)


# Regression 4: pro-poor policies and PT pro-poor policies
reg4i <- lm(ef$ideo_imp~ef$ideo_imp_lag+ef$var_pob+ef$pt+I(ef$var_pob*ef$pt)+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg4i)


# Regression 5: pooled regression
reg5i <- lm(ef$ideo_imp~ef$ideo_imp_lag+ef$closeness+ef$ideo_pres+ef$close_ideopres+ef$supp_pres+I(ef$supp_pres*ef$ideo_pres)+ef$coal_pres+ef$IDHM+ef$fragf+ef$comp+ef$turnout+ef$pol_pi+ef$var_pob+ef$ln_voter+ef$national+ef$year, data = ts(ef, start = c(1994,1), end = c(2018,13), frequency = 13))
summary(reg5i)




######################################################################################################
######################################################################################################
######################################################################################################
########################### GLS Regressions with BLS D.V.  ###########################################
######################################################################################################
######################################################################################################
######################################################################################################


library(nlme)

# Regression 1: imcumbent spillover effect
reg1.gls <- gls(ideo_na~ideo_na_lag+closeness+close_ideopres+ideo_pres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+ln_voter+national+year, data = ez , na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg1.gls)


# Regression 2: socioeconomic modernization
reg2.gls <- gls(ideo_na~ideo_na_lag+IDHM+ln_voter+national+year, data = ez, na.action="na.omit",  method="ML", correlation = corAR1(form = ~1|year))
summary(reg2.gls)


# Regression 3: political modernization (pluralism)
reg3.gls <- gls(ideo_na~ideo_na_lag+fragf+comp+turnout+pol_pi+ln_voter+national+year, data = ez , na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg3.gls)


# Regression 4: pro-poor policies and PT pro-poor policies
reg4.gls <- gls(ideo_na~ideo_na_lag+var_pob+pt+I(var_pob*pt)+ln_voter+national+year, data = ez, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg4.gls)


# Regression 5: pooled regression
reg5.gls <- gls(ideo_na~ideo_na_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ez, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg5.gls)


######################################################################################################
######################################################################################################
######################################################################################################
########################### GLS Regressions with IMPUTED D.V. ########################################
######################################################################################################
######################################################################################################
######################################################################################################


# Regression 1: imcumbent spillover effect
reg1i.gls <- gls(ideo_imp~ideo_imp_lag+closeness+close_ideopres+ideo_pres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+ln_voter+national+year, data = ez , na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg1i.gls)


# Regression 2: socioeconomic modernization
reg2i.gls <- gls(ideo_imp~ideo_imp_lag+IDHM+ln_voter+national+year, data = ez, na.action="na.omit",  method="ML", correlation = corAR1(form = ~1|year))
summary(reg2i.gls)


# Regression 3: political modernization (pluralism)
reg3i.gls <- gls(ideo_imp~ideo_imp_lag+fragf+comp+turnout+pol_pi+ln_voter+national+year, data = ez , na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg3i.gls)


# Regression 4: pro-poor policies and PT pro-poor policies
reg4i.gls <- gls(ideo_imp~ideo_imp_lag+var_pob+pt+I(var_pob*pt)+ln_voter+national+year, data = ez, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg4i.gls)


# Regression 5: pooled regression
reg5i.gls <- gls(ideo_imp~ideo_imp_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ez, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(reg5i.gls)



#########################################################################################################
#########################################################################################################
#########################################################################################################
##################################### Fig 5. PLOTS of the COEFFICIENTS ##################################
#########################################################################################################
#########################################################################################################
#########################################################################################################


library(dotwhisker)
library(broom)
library(dplyr)
library(coefplot)
library(lmtest)
library(arm)


vnm <- c("Ideo. Prox. Mayor/Pres.","Presid. Ideology","Proximity * Pres. Ideology","Pres. Approval","Approval * Ideology Pres.", "Political Alignment", "Human Development Index","Party Fragmentation","Electoral Competition","Turnout","Polarization","Variation in Poverty Rate","Size of Electorate","Concurrent Elections","Year")

model <- reg5i.gls
varnm <- names(coef(model))
varnm <- varnm[3:length(varnm)]
slopes <- coef(summary(model))[varnm, 1]  #' slopes
ses <- coef(summary(model))[varnm, 2]  #' SEs
ci <- confint(model, varnm, level = 0.99)


par(mar=c(4,11,1,1))
plot(NA, xlim = c(-0.25, 0.5), ylim = c(0, length(varnm)), xlab = "Impact on the D.V. (Z-Score)", ylab = "", yaxt = "n", bty="n")

axis(2, 1:length(varnm), vnm, las = 2)

abline(v = 0, col = "gray")

points(slopes, 1:length(slopes), pch = 19, col = "black")


segments(ci[,1], 1:length(slopes), ci[,2], 1:length(slopes), col = "black", lwd = 1)




#######################################################################################################################
#######################################################################################################################
#######################################################################################################################
######################################### Robustness check with LAGGED DV #############################################
#######################################################################################################################
#######################################################################################################################
#######################################################################################################################

res <- read.csv("Data/Prop_imputed.csv", header=T, sep=";", dec=",")



res$x <- res$x*100

res$cods <- 1
res$cods[res$x>0 & res$x<=10] <- 2
res$cods[res$x>10 & res$x<=20] <- 3
res$cods[res$x>20 & res$x<=30] <- 4
res$cods[res$x>30 & res$x<=40] <- 5
res$cods[res$x>40] <- 6

ta <- table(res$cods, res$year)
ta <- rbind(ta,colSums(ta))


for (i in 1:12) ta[1:6,i] <- round(ta[1:6,i]/ta[7,i]*100,1)


ee <- merge(ef,res, by=c("cod_tse", "year"), all.x=T)

rm(ef,res,ta,ez)

for(i in 1:20) gc()

library(nlme) 

# Regression 5: pooled regression

# No imputation
cods <- 1
ea <- ee[ee$cods%in%cods,]
pooled0 <- gls(ideo_imp~ideo_imp_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ea, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(pooled0)

# until 10% imputation
cods <- c(1,2)
ea <- ee[ee$cods%in%cods,]
pooled10 <- gls(ideo_imp~ideo_imp_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ea, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(pooled10)


# until 20% imputation
cods <- c(1,2,3)
ea <- ee[ee$cods%in%cods,]
pooled20 <- gls(ideo_imp~ideo_imp_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ea, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(pooled20)

# until 30% imputation
cods <- c(1,2,3,4)
ea <- ee[ee$cods%in%cods,]
pooled30 <- gls(ideo_imp~ideo_imp_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ea, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(pooled30)

# until 40% imputation
cods <- c(1,2,3,4,5)
ea <- ee[ee$cods%in%cods,]
pooled40 <- gls(ideo_imp~ideo_imp_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ea, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(pooled40)

# more than 40% imputation
cods <- c(1,2,3,4,5,6)
pooledm40 <- gls(ideo_imp~ideo_imp_lag+closeness+ideo_pres+close_ideopres+supp_pres+I(supp_pres*ideo_pres)+coal_pres+IDHM+fragf+comp+turnout+pol_pi+var_pob+ln_voter+national+year, data = ee, na.action="na.omit",  method="ML", correlation = corAR1(form = ~ 1|year))
summary(pooledm40)


rm(list=ls())




